function y=multireg(Y,X,v,V,B_bar,A)

%This is a function to draw sample from a posterior of a multivariate
%regression model with prior
% Sigma ~ IW(v,V)
% Beta = vec(B)| Sigma ~ N(vec(B_bar), kron(Sigma, inv(A)))

%Y is n*m, n is the number of observation, m is the number of equations
%X is n*k, k is the number of parameters in each equation
%B is k*m

n=size(Y,1);
m=size(Y,2);
k=size(X,2);


%Step 1: draw Sigma 
[Q_A R_A]=qr(A,0);
W=[X;R_A];
Z=[Y;R_A*B_bar];
R_ww=chol(W'*W);
IR_ww=inv(R_ww);

Beta=IR_ww*IR_ww'*W'*Z;
E=Z-W*Beta;
S=E'*E;

sigma=iwhish(V+S,v+n);


%Step 2: draw beta conditional on sigma
C=chol(inv(sigma));
CI=inv(C);

z=randn(m*k,1);

beta=Beta+IR_ww*reshape(z,k,m)*CI';

y=[beta;sigma];